---
title: "01_figures"
date: "9/28/2021"
output: html_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(srvyr)
library(haven)
library(questionr)

## This code is ordered as follow:
## 1. Density plots: Figures 5, 9-13
## 2. Barplot: Figures 1-4, 6-7

```


###############################################################################################
###############################################################################################
# Code for density plots

## read dta
```{r}
cep<- read_dta("cep84.dta")

cep_svy <- cep %>% 
  as_survey_design(strata = strata_cep, id = UMP, weights = PONDERADOR, nest = T)

```

## Figures 

```{r}
cep <- cep %>% mutate(lib_orden=ifelse(ESP_4<20, ESP_4, NA), 
                                  tipo2=6-tipo,  ig_esf=ifelse(ESP_11<20, ESP_11, NA), 
                      resp_sust=ifelse(ESP_12<20, ESP_12, NA), 
                      tipo_s=factor(tipo2, levels=c(1, 2, 3, 4, 5),
                             labels=c("Undecided", "Opponent",
                                      "Supporter", "Weak protester", "Strong protester" ), ordered = TRUE),
                      tipo_s2=factor(tipo, levels=c(1, 2, 3, 4, 5),
                                    labels=c("Strong protester", "Weak protester",  "Supporter", "Opponent",
                                             "Undecided"), ordered = TRUE))

```

### Figure 5

#### Figure 5_a equality vs effort

```{r}

# equality vs effort

g5_a <- ggplot(cep, aes(ig_esf, colour=tipo_s2,linetype=tipo_s2, shape=tipo_s2, weight=PONDERADOR)) +
  geom_density(alpha = 0.3) + theme_classic()  + scale_x_continuous( limits = c(0,10)) + 
  xlab("") + ylab("Density") +
scale_y_continuous(expand= c(0,0), limits = c(0, 0.3)) + 
  scale_color_manual(values=c("grey0","grey20","gray30", "gray60", "grey50"))+
    theme(axis.text.x=element_text(size=12), legend.position=c(.9,.8),
           plot.title = element_text(lineheight=1, face="bold", hjust = .5),
        plot.caption=element_text(hjust = 0.5),
        legend.text = element_text(size = 12)) +
  labs(colour  = "", linetype = "", shape = "") +
scale_x_continuous(breaks = 1:10,
    labels = c("1 \n Equality", "2", "3", "4", "5", "6", "7", "8", "9", "10 \n  Individual \n effort"))

ggsave("./figures/Figure_5a.pdf", g5_a, width = 9, height = 5)

```

#### fig 5b 


```{r}


g5_b <- ggplot(cep, aes(resp_sust, colour=tipo_s2,linetype=tipo_s2, shape=tipo_s2, weight=PONDERADOR)) +
  geom_density(alpha = 0.3) + theme_classic()  +  
  xlab("") + ylab("Density") +
scale_y_continuous(expand= c(0,0), limits = c(0, 0.3)) + 
  scale_color_manual(values=c("grey0","grey20","gray30", "gray60", "grey50"))+
    theme(axis.text.x=element_text(size=12), legend.position=c(.9,.8),
           plot.title = element_text(lineheight=1, face="bold", hjust = .5),
        plot.caption=element_text(hjust = 0.5, size=12),
        legend.text = element_text(size = 12)) +
  labs(colour  = "", linetype = "", shape = "") +
scale_x_continuous(breaks = 1:10,
    labels = c("1 \n Government", "2", "3", "4", "5", "6", "7", "8", "9", "10 \n  People \n themselves")) 

ggsave("./figures/Figure_5b.pdf", g5_b, width = 9, height = 5)
```


### Figure 9
```{r}
g9 <- ggplot(cep, aes(lib_orden, colour=tipo_s2, linetype=tipo_s2, shape=tipo_s2, weight=PONDERADOR)) +
  geom_density(alpha = 0.3) + theme_classic()  +
  xlab("") + ylab("Density") +
scale_y_continuous(expand= c(0,0), limits = c(0, 0.3)) + labs(colour = "")   +
  scale_color_manual(values=c(
                              "grey0","grey20","gray30", "gray60", "grey50"))+
    theme(axis.text.x=element_text(size=12), legend.position=c(.9,.8),
           plot.title = element_text(lineheight=1, face="bold", hjust = .5),
        plot.caption=element_text(hjust = 0.5),
        legend.text = element_text(size = 12)) +
   
  labs(colour  = "", linetype = "", shape = "") +
scale_x_continuous(breaks = 1:10,
    labels = c("1 \n Public and  \n private liberties", "2", "3", "4", "5", "6", "7", "8", "9", "10 \n  Public \n order")) 

ggsave("./figures/Figure_9.pdf", g9, width = 8.5, height = 5)

```


### Figure 10
```{r}


res_dem<- cep_svy %>% mutate(e1=as.character(e1), tipo=as.character(tipo)) %>%  group_by(tipo,e1) %>% 
  summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit"))  %>% filter(e1=="1") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=1)


res_fun<- cep_svy %>% mutate(e2=as.character(e2), tipo=as.character(tipo)) %>%  group_by(tipo,e2) %>% 
  summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit"))  %>% filter(e2=="3") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=2)

resbdem<- rbind(res_dem,  res_fun)  %>% 
          mutate(prop_low=ifelse(prop_low<0, 0, prop_low), 
                       prop=100*prop, prop_low=100*prop_low, prop_upp=100*prop_upp,
                       lege=factor(leg2, levels=c(1, 2),
                            labels=c("Democracy is preferable to any other kind of government", 
                                      "Negative evaluation of how democracy is working"), ordered = TRUE))

dodge <- position_dodge(width=0.9)

g10 <- ggplot(resbdem) + geom_bar(aes(x=tipo, y=prop, fill=lege),width = .8, colour="black", stat="identity", position = dodge) +
  geom_linerange(aes(x=tipo, ymin=prop_low, ymax=prop_upp, fill=lege), stat="identity", position = dodge, width = 0.25) +
  scale_x_discrete(name ="",  labels=c("1"="Strong protester","2"= "Weak protester", "3"="Supporter", "4"="Opponent", "5"="Undecided")) +
  scale_fill_manual(values=c("grey80","grey20")) + scale_y_continuous(expand=c(0,0), limits=c(0,100), name="Percent", breaks=seq(0,100, 10))+
  guides(fill=guide_legend(title=NULL)) + theme_classic() + 
  theme(axis.text.x=element_text(size=10), legend.position=c(.68,.9),
        plot.caption=element_text(hjust = 0.5), axis.line.x = element_blank(),
        plot.title = element_text(lineheight=1, face="bold", hjust = .5),
         plot.subtitle = element_text(hjust = 0.5, size=8)) + 
  labs(caption  = "*Black bars show 95% confidence intervals")
       #, title="Perceptions of democracy")

ggsave("./figures/Figure_10.pdf", g10, width = 8.5, height = 5)
```



### Figure 11
```{r}

##4.	Los partidos políticos
##12.	El Gobierno
##13.	El Congreso




res_conf1<- cep_svy %>% mutate( tipo=as.character(tipo)) %>%  group_by(tipo,conf12) %>% 
  summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit"))  %>% filter(conf12=="1") %>% 
  select(tipo, prop, prop_low, prop_upp) %>% mutate( conf=1)

res_conf2<- cep_svy %>% mutate( tipo=as.character(tipo)) %>%  group_by(tipo,conf13) %>% 
  summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit"))  %>% filter(conf13=="1") %>% 
  select(tipo, prop, prop_low, prop_upp) %>% mutate( conf=2)


res_conf3<- cep_svy %>% filter(tipo!=2) %>% mutate( tipo=as.character(tipo)) %>%  group_by(tipo,conf4) %>% 
  summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit"))  %>% filter(conf4=="1") %>% 
  select(tipo, prop, prop_low, prop_upp) %>% mutate( conf=3)

tipo2 <- tibble(tipo="2", prop=0, prop_low=0, prop_up=0, conf=3)

res_conf3 <- rbind(res_conf3, tipo2)


res_conf<- rbind(res_conf1, res_conf2, res_conf3)  %>% 
  mutate(prop_low=ifelse(prop_low<0, 0, prop_low), prop=100*prop, 
         prop_low=100*prop_low, prop_upp=100*prop_upp,
         lege=factor(conf, levels=c(1, 2, 3),
          labels=c("Government", "Congress", "Political parties"), ordered = TRUE))


dodge <- position_dodge(width=0.9)
g11 <- ggplot(res_conf) + geom_bar(aes(x=tipo, y=prop, fill=lege),width = .8, colour="black", stat="identity", position = dodge) +
  geom_linerange(aes(x=tipo, ymin=prop_low, ymax=prop_upp, fill=lege), stat="identity", position = dodge, width = 0.25) +
  scale_x_discrete(name ="",  labels=c("1"="Strong protester","2"= "Weak protester", "3"="Supporter", "4"="Opponent", "5"="Undecided")) +
  scale_fill_manual(values=c("grey0","grey40", "grey100")) + scale_y_continuous(expand=c(0,0), limits=c(0,20), name="Percent", breaks=seq(0,20, 5))+
  guides(fill=guide_legend(title=NULL)) + theme_classic() + 
  theme(axis.text.x=element_text(size=10), legend.position=c(.85,.85),
        plot.caption=element_text(hjust = 0.5), axis.line.x = element_blank(),
        plot.title = element_text(lineheight=1, face="bold", hjust = .5),
         plot.subtitle = element_text(hjust = 0.5, size=8)) + 
  labs(caption  = "*Black bars show 95% confidence intervals",
       title="")

ggsave("./figures/Figure_11.pdf", g11, width = 8.5, height = 5)
```





###  Fig 12: JUSTIFIC VIOLENCIA
  
```{r}


  res_barric <- cep_svy %>% mutate(e42=as.character(ifelse(e41_barric==1, 1, 0)), 
                                   tipo=as.character(tipo)) %>% group_by(tipo, e42) %>% 
    summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit")) %>% 
  filter(e42=="1") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=2)
  
  res_saq <- cep_svy %>% mutate(e42=as.character(ifelse(e41_saqueo==1, 1, 0)), tipo=as.character(tipo)) %>% 
    group_by(tipo, e42) %>% 
    summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit")) %>% filter(e42=="1") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=3)
  
  
  res_evad<- cep_svy %>% mutate(e42=as.character(ifelse(e41_evad==1, 1, 0)), tipo=as.character(tipo)) %>% 
    group_by(tipo, e42) %>% 
    summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit"))  %>% 
    filter(e42=="1") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=1)
  
  
  res_tear <- cep_svy %>% mutate(e42=as.character(ifelse(e42_1==1, 1, 0)), tipo=as.character(tipo)) %>% 
    group_by(tipo, e42) %>% 
    summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit")) %>% 
    filter(e42=="1") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=4)
  
  res_f<- cep_svy %>% mutate(e42=as.character(ifelse(e42_2==1, 1, 0)), tipo=as.character(tipo)) %>% 
    group_by(tipo, e42) %>% 
    summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit")) %>% filter(e42=="1") %>% 
    select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=5)
  
  
  resbase<- rbind(res_evad, res_barric, res_saq , res_tear, res_f)  %>% 
    mutate(prop_low=ifelse(prop_low<0, 0, prop_low),prop=100*prop, prop_low=100*prop_low, prop_upp=100*prop_upp,
            lege=factor(leg2, levels=c(1, 2, 3, 4, 5),
                        labels=c("Evade transport fee to protest", 
                                "Participate in barricades to protest",
                                  "Participate in looting to protest",
                                "Police use force against a violent protester",
                                 "Police use tear gas to control violence at rally"), ordered = TRUE))
  
  dodge <- position_dodge(width=0.9)
  
 
   g12 <-    ggplot(resbase) +
    geom_bar(aes(x=tipo, y=prop, fill=lege), width = .8, colour="black", stat="identity", position = dodge) +
    geom_linerange(aes(x=tipo, ymin=prop_low, ymax=prop_upp, fill=lege), stat="identity", position = dodge, width = 0.25) +
    scale_x_discrete(name ="",  labels=c("1"="Strong protester","2"= "Weak protester", "3"="Supporter", "4"="Opponent", "5"="Undecided")) +
    scale_fill_manual(values=c("grey0","grey25","gray50", "gray80", "grey100")) + 
    scale_y_continuous(expand=c(0,0), limits=c(0,100), name="Percent", breaks=seq(0,100, 10))+
    guides(fill=guide_legend(title=NULL)) + theme_classic() + 
    theme(axis.text.x=element_text(size=12), legend.position=c(.8,.80),
            plot.title = element_text(lineheight=1, face="bold", hjust = .5),
            plot.caption=element_text(hjust = 0.5),
            plot.subtitle = element_text(hjust = 0.5, size=8),
            axis.line.x = element_blank(),
            legend.text = element_text(size=10)) + 
    labs(      caption  = "*Black bars show 95% confidence intervals")
   
    ggsave("./figures/Figure_12.pdf", g12, width = 8.5, height = 5)
    
```


###  Fig 13: abuse human rights caused fire : police & military

```{r}
cep84_resp <- read_dta("cep84_resp.dta")
cep2_svy <- cep84_resp %>% 
  as_survey_design(strata = strata_cep, id = UMP, weights = PONDERADOR, nest = T)

rescar <- cep_svy %>% mutate(e43=as.character(ifelse(e43==1 | e44==1, 1, 0)), tipo=as.character(tipo)) %>% group_by(tipo, e43) %>% 
  summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit")) %>% filter(e43=="1") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=1)

resfire <- cep2_svy %>% mutate(tipo=as.character(tipo)) %>% group_by(tipo, d39carab) %>% 
  summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit"))  %>% filter(d39carab=="1") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=2)

ressaq <- cep2_svy %>% mutate(tipo=as.character(tipo)) %>% group_by(tipo, d40carab) %>% 
  summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit"))  %>% filter(d40carab=="1") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=3)


res_tc<- cep_svy %>% mutate( tipo=as.character(tipo)) %>%  group_by(tipo,conf11) %>% 
  summarize(prop=survey_mean(na.rm = T, vartype = "ci", prop_method = "logit"))  %>% filter(conf11=="1") %>% select(tipo, prop, prop_low, prop_upp) %>% mutate( leg2=4)



resbase<- rbind(rescar,  resfire, ressaq, res_tc)  %>% 
    mutate(prop_low=ifelse(prop_low<0, 0, prop_low), 
           prop=100*prop, prop_low=100*prop_low, prop_upp=100*prop_upp,
           lege=factor(leg2, levels=c(1, 2, 3, 4),
           labels=c("Police or military abused human rights frequently or very frequently", 
                    "Believes police or military are responsible for the fires",
                    "Believes police or military are responsible for looting",
                     "Trust the police"), ordered = TRUE))


dodge <- position_dodge(width=.9)
g13<- ggplot(resbase) + geom_bar(aes(x=tipo, y=prop, fill=lege), width = .8, colour="black", stat="identity", position = dodge) +
  geom_linerange(aes(x=tipo, ymin=prop_low, ymax=prop_upp, fill=lege), stat="identity", position = dodge, width = 0.25) +
  scale_x_discrete(name ="",  labels=c("1"="Strong protester","2"= "Weak protester", "3"="Supporter", "4"="Opponent", "5"="Undecided")) +
  scale_fill_manual(values=c("grey0","grey25","grey55", "grey100")) + scale_y_continuous(expand=c(0,0), limits=c(0,100), name="Percent", breaks=seq(0,100, 10))+
guides(fill=guide_legend(title=NULL)) + theme_classic() + 
  theme(axis.text.x=element_text(size=10), legend.position=c(.72,.85),
         plot.caption=element_text(hjust = 0.5),
        plot.title = element_text(lineheight=1, face="bold", hjust = .5)) + 
  labs(title="", caption  = "*Black bars show 95% confidence intervals") 

ggsave("./figures/Figure_13.pdf", g13, width = 8.5, height = 5)
```




###############################################################################################
###############################################################################################
# Code for barplots


```{r}

d <- cep


## Creating variable "tipo" with our typology
d$tipo<-as.factor(d$tipo)
levels(d$tipo)<-c("Strong protester", "Weak protester", "Supporter", "Opponent", "Undecided")


## sample size by type (weighted by survey weights: d$PONDERADOR)
n_types<-wtd.table(d$tipo, weights=d$PONDERADOR)

```



### Figure 1: Distribution of groups
```{r}
## proportion by type (weighted by survey weights: PONDERADOR)
# n=1496
p<-wtd.table(d$tipo, weights=d$PONDERADOR)/1496 


pdf(file = "./figures/fig_1.pdf",   # The directory you want to save the file in
    width = 8.5, # The width of the plot in inches
    height = 4)
barplot(100*p, col="gray50",names.arg=c("Strong protester", "Weak protester", "Supporter", "Opponent", "Undecided"),
        ylab = "Percent",ylim=c(0,55),xlim=c(0,6))
std.error<-100*sqrt(p*(1-p)/1496)
mtext("*Black bars show 95% confidence intervals",cex=.8,side=1,line=3.7, las=0)
segments(seq(.65,by=1.2,length.out =5),100*p-1.96*std.error,seq(.65,by=1.2,length.out =5),100*p+1.96*std.error,lwd=2)
text(seq(from=.65,to=6,by=1.2), 100*p+6, paste(round(100*p),"%", sep=""))
dev.off()


```

### Figure 2: Age distribution by group
```{r}

## Create variables for young and old
d$young<-ifelse(d$EDAD_T<3,1,0)
d$old<-ifelse(d$EDAD_T==5,1,0)
#means by groups
tapply(d$young,d$tipo, mean, weights.arg=d$PONDERADOR)
tapply(d$old,d$tipo, mean, weights.arg=d$PONDERADOR)

## create weighted table by age group and type
age_type<-rbind(tapply(seq(along=d$young), list(d$tipo), function(i, x=d$young, w=d$PONDERADOR)
  weighted.mean(x[i], w[i])) , ##distribution of young by type, weighted
  tapply(seq(along=d$old), list(d$tipo), function(i, x=d$old, w=d$PONDERADOR)
    weighted.mean(x[i], w[i]))) ##distribution of old by type, weighted

pdf(file = "./figures/fig_2.pdf",   # The directory you want to save the file in
    width = 8.5, # The width of the plot in inches
    height = 4)
barplot(100*age_type,beside=T,col=c("gray20","gray80"),ylim=c(0,70),
        names.arg=c("Strong protester", "Weak protester", "Supporter", "Opponent", "Undecided"),ylab="Percent")
legend(col=c("gray20","gray80"),pch=15, legend=c("Young (18-34)", "Old (55 or more)"),"topright", cex=1)

## Standard errors: sqrt(p*(1-p))/n_i
st.error<-100*sqrt((age_type*(1-age_type))/rbind(n_types,n_types))
segments(c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),100*age_type-1.96*st.error,c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),100*age_type+1.96*st.error)
mtext("*Black bars show 95% confidence intervals",cex=.8,side=1,line=3.7, las=0)
dev.off()


```


### Figure 3: Education by group

```{r}
#means by groups
tapply(d$higher_ed,d$tipo, mean)
tapply(d$no_high_school,d$tipo, mean, weights.arg=d$PONDERADOR)

educ_type<-rbind(tapply(seq(along=d$higher_ed), list(d$tipo), function(i, x=d$higher_ed, w=d$PONDERADOR)
  weighted.mean(x[i], w[i])) ,
  tapply(seq(along=d$no_high_school), list(d$tipo), function(i, x=d$no_high_school, w=d$PONDERADOR)
    weighted.mean(x[i], w[i])))

pdf(file = "./figures/fig_3.pdf",   # The directory you want to save the file in
    width = 8.5, # The width of the plot in inches
    height = 4)
barplot(100*educ_type,beside=T,col=c("gray20","gray80"),ylim=c(0,60),
        names.arg=c("Strong protester", "Weak protester", "Supporter", "Opponent", "Undecided"),ylab="Percent")
legend(col=c("gray20","gray80"),pch=15, legend=c("With higher education", "Without complete high school"),"topleft", cex=1)
st.error<-100*sqrt((educ_type*(1-educ_type))/rbind(n_types, n_types))
segments(c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),100*educ_type-1.96*st.error,c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),100*educ_type+1.96*st.error)
mtext("*Black bars show 95% confidence intervals",cex=.8,side=1,line=3.7, las=0)
dev.off()


```


### Figure 4: political position by group
```{r}
##Create variable for political position from position in 1-10 ideological axis
d$pospol<-ifelse(d$DS_P20<5,1,ifelse(d$DS_P20<7,2,ifelse(d$DS_P20<11,3,4)))
d$pospol<-as.factor((d$pospol))
levels(d$pospol)<-c("Left", "Center", "Right", "No ID")

pospol_type<-wtd.table(d$pospol,d$tipo,weights=d$PONDERADOR)/rbind(n_types,n_types,n_types,n_types)

pdf(file = "./figures/fig_4.pdf",   # The directory you want to save the file in
    width = 8.5, # The width of the plot in inches
    height = 4)
barplot(100*pospol_type, beside=T,ylim=c(0,80),col=c("gray0","gray25","gray50", "gray100"),
        ylab="Percent")
st.error<-100*sqrt((pospol_type*(1-pospol_type))/rbind(n_types,n_types,n_types,n_types))
segments(c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5, 14.5,16.5,17.5,18.5,19.5,21.5,22.5,23.5,24.5),
         100*pospol_type-1.96*st.error,c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5, 14.5,16.5,17.5,18.5,19.5,21.5,22.5,23.5,24.5),100*pospol_type+1.96*st.error)
mtext("*Black bars show 95% confidence intervals",cex=.8,side=1,line=3.7, las=0)
legend(col=c("gray0","gray25","gray50", "gray100"),legend =c("Left", "Centre", "Right", "No ID"),pch=15,"top",horiz=T )
dev.off()



```


### Figure 6: Salary Considered Appropriate vs. Perceived Real Salary for Different Occupations by Group

```{r}
## Means and medians by group
tapply(d$ESP_14A1[d$ESP_14A1!=99],d$tipo[d$ESP_14A1!=99],median,weights.arg=d$PONDERADOR[d$ESP_14A1!=99])/1000
tapply(d$ESP_14B1[d$ESP_14B1!=99],d$tipo[d$ESP_14B1!=99],median,weights.arg=d$PONDERADOR[d$ESP_14B1!=99])/1000
tapply(d$ESP_14C1[d$ESP_14C1!=99],d$tipo[d$ESP_14C1!=99],median,weights.arg=d$PONDERADOR[d$ESP_14C1!=99])/1000

tapply(d$ESP_14A1[d$ESP_14A1!=99],d$tipo[d$ESP_14A1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_14A1!=99])/1000
tapply(d$ESP_14B1[d$ESP_14B1!=99],d$tipo[d$ESP_14B1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_14B1!=99])/1000
tapply(d$ESP_14C1[d$ESP_14C1!=99],d$tipo[d$ESP_14C1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_14C1!=99])/1000

tapply(d$ESP_15A1[d$ESP_15A1!=99],d$tipo[d$ESP_15A1!=99],median,weights.arg=d$PONDERADOR[d$ESP_15A1!=99])/1000
tapply(d$ESP_15B1[d$ESP_15B1!=99],d$tipo[d$ESP_15B1!=99],median,weights.arg=d$PONDERADOR[d$ESP_15B1!=99])/1000
tapply(d$ESP_15C1[d$ESP_15C1!=99],d$tipo[d$ESP_15C1!=99],median,weights.arg=d$PONDERADOR[d$ESP_15C1!=99])/1000

tapply(d$ESP_15A1[d$ESP_15A1!=99],d$tipo[d$ESP_15A1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_15A1!=99])/1000
tapply(d$ESP_15B1[d$ESP_15B1!=99],d$tipo[d$ESP_15B1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_15B1!=99])/1000
tapply(d$ESP_15C1[d$ESP_15C1!=99],d$tipo[d$ESP_15C1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_15C1!=99])/1000



## Function for weighted variance 
weighted.var <- function(x, w, na.rm = FALSE) {
  if (na.rm) {
    w <- w[i <- !is.na(x)]
    x <- x[i]
  }
  sum.w <- sum(w)
  sum.w2 <- sum(w^2)
  mean.w <- sum(x * w) / sum(w)
  (sum.w / (sum.w^2 - sum.w2)) * sum(w * (x - mean.w)^2, na.rm =
                                       na.rm)
}


```


#### Figures 6a-c
```{r}

## weighted means and SE excluding "don't know" (99)
worker<-rbind(tapply(d$ESP_14C1[d$ESP_14C1!=99],d$tipo[d$ESP_14C1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_14C1!=99])/1000,
              tapply(d$ESP_15C1[d$ESP_15C1!=99],d$tipo[d$ESP_15C1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_15C1!=99])/1000)

worker_se<-rbind(sqrt(tapply(d$ESP_14C1[d$ESP_14C1!=99],d$tipo[d$ESP_14C1!=99],weighted.var,w=d$PONDERADOR[d$ESP_14C1!=99] )/n_types)/1000,
                 sqrt(tapply(d$ESP_15C1[d$ESP_15C1!=99],d$tipo[d$ESP_15C1!=99],weighted.var,w=d$PONDERADOR[d$ESP_15C1!=99] )/n_types)/1000)


pdf(file = "./figures/fig_6a.pdf",   # The directory you want to save the file in
    width = 8.5, # The width of the plot in inches
    height = 4)
barplot(worker, beside=T, ylim=c(0,2500),ylab="Mean in thousands of pesos",col=c("gray20","gray80"))
legend(col=c("gray20","gray80"),pch=15, legend=c("Earns","Should earn"),"topright",horiz=T)
segments(c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),worker-1.96*worker_se,c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),worker+1.96*worker_se)
mtext("*Black bars show 95% confidence intervals",cex=.8,side=1,line=3.7, las=0)
dev.off()


#Figure 6B: doctor
doctor<-rbind(tapply(d$ESP_14A1[d$ESP_14A1!=99],d$tipo[d$ESP_14A1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_14A1!=99])/1000,
              tapply(d$ESP_15A1[d$ESP_15A1!=99],d$tipo[d$ESP_15A1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_15A1!=99])/1000)

doctor_se<-rbind(sqrt(tapply(d$ESP_14A1[d$ESP_14A1!=99],d$tipo[d$ESP_14A1!=99],weighted.var,w=d$PONDERADOR[d$ESP_14A1!=99] )/n_types)/1000,
                 sqrt(tapply(d$ESP_15A1[d$ESP_15A1!=99],d$tipo[d$ESP_15A1!=99],weighted.var,w=d$PONDERADOR[d$ESP_15A1!=99] )/n_types)/1000)

pdf(file = "./figures/fig_6b.pdf",   # The directory you want to save the file in
    width = 8.5, # The width of the plot in inches
    height = 4)
barplot(doctor, beside=T, ylim=c(0,25000),ylab="Mean in thousands of pesos",col=c("gray20","gray80"))
legend(col=c("gray20","gray80"),pch=15, legend=c("Earns","Should earn"),"topright",horiz=T)
segments(c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),doctor-1.96*doctor_se,c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),doctor+1.96*doctor_se)
mtext("*Black bars show 95% confidence intervals. Note that the scale is 10x than in Figure 6A.",cex=.8,side=1,line=3.7, las=0)
dev.off()

#Figure 6C: CEO
CEO<-rbind(tapply(d$ESP_14B1[d$ESP_14B1!=99],d$tipo[d$ESP_14B1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_14B1!=99])/1000,
           tapply(d$ESP_15B1[d$ESP_15B1!=99],d$tipo[d$ESP_15B1!=99],mean,weights.arg=d$PONDERADOR[d$ESP_15B1!=99])/1000)

CEO_se<-rbind(sqrt(tapply(d$ESP_14B1[d$ESP_14B1!=99],d$tipo[d$ESP_14B1!=99],weighted.var,w=d$PONDERADOR[d$ESP_14B1!=99] )/n_types)/1000,
              sqrt(tapply(d$ESP_15B1[d$ESP_15B1!=99],d$tipo[d$ESP_15B1!=99],weighted.var,w=d$PONDERADOR[d$ESP_15B1!=99] )/n_types)/1000)

pdf(file = "./figures/fig_6c.pdf",   # The directory you want to save the file in
    width = 8.5, # The width of the plot in inches
    height = 4)
barplot(CEO, beside=T, ylim=c(0,25000),ylab="Mean in in thousands of pesos",col=c("gray20","gray80"))
legend(col=c("gray20","gray80"),pch=15, legend=c("Earns","Should earn"),"topright",horiz=T)
segments(c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),CEO-1.96*CEO_se,c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5),CEO+1.96*CEO_se)
mtext("*Black bars show 95% confidence intervals. Note that the scale is 10x than in Figure 6A.",cex=.8,side=1,line=3.7, las=0)
dev.off()

```




### Figure 7: Two Main Causes of Economic Success by Group (per cent)


```{r}
causas_riqueza<-wtd.table(as.factor(d$ESP_10_1),d$tipo,weights=d$PONDERADOR)+wtd.table(as.factor(d$ESP_10_2),d$tipo,weights=d$PONDERADOR)
causas_riqueza<-prop.table(causas_riqueza,2) ##column percentage
```

```{r}


## Sum up first and second main reasons


## order by most mentioned
order(prop.table((wtd.table(as.factor(d$ESP_10_1),weights=d$PONDERADOR)+wtd.table(as.factor(d$ESP_10_2),weights=d$PONDERADOR))))
## picking the 5 most mentioned
causas_riqueza<-causas_riqueza[c(6,4,5,1,7),]
causas_riqueza<-t(causas_riqueza) # transpose
st.error<-100*sqrt((causas_riqueza*(1-causas_riqueza))/cbind(n_types,n_types,n_types,n_types,n_types))

# Responses short names
colnames(causas_riqueza)<-c("Education", "Responsible work","Networks","Personal initiative", "Parents' situation")


pdf(file = "./figures/fig_7.pdf",   # The directory you want to save the file in
    width = 8.5, # The width of the plot in inches
    height = 4)
barplot(100*causas_riqueza, beside=T,col=c("gray0","gray25","gray50", "gray75","gray100"),
        ylab="Percent",ylim=c(0,40), cex.names = .8)
legend(col=c("gray0","gray25","gray50", "gray75","gray100"),legend =c("Strong protester","Weak protester","Supporter","Opponent", "Undecided"),pch=15,"topright", cex = .8 )
segments(c(1.5,2.5,3.5,4.5,5.5,7.5,8.5,9.5,10.5,11.5,13.5,14.5, 15.5,16.5,17.5,19.5,20.5,21.5,22.5,23.5,25.5,26.5,27.5,28.5,29.5,31.5,32.5,33.5,34.5,35.5),
         100*causas_riqueza-1.96*st.error,
         c(1.5,2.5,3.5,4.5,5.5,7.5,8.5,9.5,10.5,11.5,13.5,14.5, 15.5,16.5,17.5,19.5,20.5,21.5,22.5,23.5,25.5,26.5,27.5,28.5,29.5,31.5,32.5,33.5,34.5,35.5),100*causas_riqueza+1.96*st.error)
mtext("*Five most frequent choices. Black bars show 95% confidence intervals",cex=.8,side=1,line=3.7, las=0)
#mtext("()",cex=.8,side=3,line=.5, las=0)
dev.off()



```


#### Causes of poverty (mentioned in footnote 7)
```{r}

order(prop.table(wtd.table(as.factor(d$ESP_9_1),weights=d$PONDERADOR)))
order(prop.table(c(wtd.table(as.factor(d$ESP_9_1),weights=d$PONDERADOR),0)+wtd.table(as.factor(d$ESP_9_2),weights=d$PONDERADOR)))
causas_pobreza<-c(wtd.table(as.factor(d$ESP_9_1),weights=d$PONDERADOR),0)+wtd.table(as.factor(d$ESP_9_2),d$tipo,weights=d$PONDERADOR)
causas_pobreza<-prop.table(causas_pobreza,2)

causas_pobreza<-causas_pobreza[c(3,8,2,10,6,4),]
causas_pobreza<-t(causas_pobreza)
st.error<-100*sqrt((causas_pobreza*(1-causas_pobreza))/cbind(n_types,n_types,n_types,n_types,n_types,n_types))

colnames(causas_pobreza)<-c("Lack of education", "Few jobs","Lazyness","Abuses", "Egoism","Lack of gov. aid")


barplot(100*causas_pobreza, beside=T,col=c("tomato","steelblue","orange1", "springgreen2","gray"),
        ylab="Percent",main="Two major causes of poverty by type",ylim=c(0,40), cex.names = .8)
legend(col=c("tomato","steelblue","orange1", "springgreen2","gray"),legend =c("Strong protester","Weak protester","Supporter","Opponent", "Undecided"),pch=15,"topright", cex = .8 )
segments(c(1.5,2.5,3.5,4.5,5.5,7.5,8.5,9.5,10.5,11.5,13.5,14.5, 15.5,16.5,17.5,19.5,20.5,21.5,22.5,23.5,25.5,26.5,27.5,28.5,29.5,31.5,32.5,33.5,34.5,35.5),
         100*causas_pobreza-1.96*st.error,
         c(1.5,2.5,3.5,4.5,5.5,7.5,8.5,9.5,10.5,11.5,13.5,14.5, 15.5,16.5,17.5,19.5,20.5,21.5,22.5,23.5,25.5,26.5,27.5,28.5,29.5,31.5,32.5,33.5,34.5,35.5),100*causas_pobreza+1.96*st.error)
mtext("*Black bars show 95% confidence intervals",cex=.8,side=1,line=3.7, las=0)
mtext("(6 most frequent choices)",cex=.8,side=3,line=.5, las=0)


```

